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We present molecular dynamics simulations of the influence of point defects on the shear elastic 
coefficients of copper. We find that vacancies do not influence these coeflicients at all, while the 
introduction of interstitials causes a large reduction in the elastic coefficients. The simulations 
establish a phase diagram of the melting temperature as a function of the density of interstitials. 
A crystal having no free surface undergoes bulk mechanical melting as a result of the vanishing of 
C = (Cii — Ci2)/2 once the specific volume reaches a critical value, equal to the experimental 
volume of liquid phase. This critical volume is history independent, in the sense that it can be 
reached either by heating the crystal or by adding defects at a constant temperature. These results 
generalize the Born model of melting for the case where point defects are present. 
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I. INTRODUCTION 

Single phase models of bulk melting usually treat this process as a destabilization of the solid due to proliferation 
of crystalline defects or excitations. Among these, Lindemann suggested that melting occurs when the amplitude of 
atomic vibrations ^exceeds a critical valueo. Born proposed that melting is connected with the vanishing of the shear 
elastic coefficientsB. Other models suggest that melting is a consequence of the proliferation of intrinsic crystalline 
defects such as vacanciesa. However, there is no clear experimental evidence supporting any one of these models, all of 
which rely on a single, type of defect or excitation. A model of melting involving interstitials coupled to shear strains 
proposed by GranatoD is an exception as it involves coupling of defects and excitations. |-| 

Recently, a dramatic decrease of the shear resistance was observed near the bcc-hcp transition of solid ''HeO. These 
observations could have been interpreted as resulting from a coupling of point defects with a phonon which softens 
near the transition. Were it not for the bcc-hcp transition, one could envisage a continuation of this process until 
the crystal would melt. Thus, we have experimental evidence of a microscopic process which can potentially lead to 
melting. However, this scenario requires a coupling of phonons and point defects, and therefore supports a picture 
similar to Granato's, in which more than one type of excitation is involved. To understand this further, it is necessary 
to investigate how point defects couple to phonons, and how this influences the shear resistance. In the following, we 
report results of molecular dynamics (MD) simulations of the influence of point defects on the elastic properties of 
copper. Copper was chosen because of the large body of experimental data to which the results of the simulations 
can be compared. We focus on the shear elastic coefficients, with the intention to test the Born model in the presence 
of point defects. 

II. SIMULATION METHOD 

We performed molecular dynamics simulations on samples containing up to 2048 atoms using the "Tight-Binding 
(TB) Potential" with parameters calculated by Ckri and RosatoLl appropriate to copper. Newton's equations were 
solved using Gear's Predictor-Corrector algorithm^. The time-step used was about 1 fSpc. For temperature control 
we used the method known as the Nose-Hoover thermostatEI. The method of BcrcndsenE£l was used to calculate the 
volume of the samples. This method does not generate trajectories in the constant NPT (here, N denotes the constant 
number of atoms, P is the constant pressure, and T is the temperature) ensemble, but it sets the mean value of the 
pressure to the desired one while predicting tjhe correct volume of the system. For ensemble average calculations we 
used the Parrinello and Rahman (PR) methodO. The statistical ensemble created by the Parrinello Rahmauj^miibined 
with Nose's thermostat was identified as the isothermal-isotension (NtT) ensemble of statistical mechanicscj'E3(here, 
t is the constant tension). We calculated the elastic constants using a fluctuation formula in both the canonical and 
in the NtT ensembles. 

In the NtT ensemble the formula for the isothermal elastic constants Cijki is 
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5{e^jeki) = kBTC-^^^JVa (2.1) 

where e is the strain and Vq is the vohime of the sample at zero strain. Calculating elastic constants using this formula 
involves two steps. First we have to find the reference matrix, h^, defining the shape of the simulation box. The 
condition which determines the shape of the box is that the sample is at zero stress. After running the simulation to 
calculate ft-o, the program should be used again to calculate the fluctuations. No t many studies have been made using 
this method because of the very slow convergence of the strain fluctuationsoEj. 

In the constant NVT ensemble (here, N denotes the constant number of atoms, V is the constant volume, and 
T is the_temperature) , MD can determine the elastic constants by using fluctuation formulas derived by Ray and 
RahmanE^I: 

V 2,N K T 

Ci-jkl = --r^<5(cr.yCrfe;) H — {SikS-ji + 6irnSjk)+ < B > (2.2) 

Kb -I Vo 

where dij is the Kronecker delta and < _B > is the ensemble average of the Born term, defined as 

B = [d^V/drl, - idV/drab)/rab)raMrabjrabkrabi (2.3) 

where ijkl are the indices of the elastic constants tensor. This formula for calculating the elastic constants is much 
less compact than the NtT formula, mainly because it contains the explicit second derivative of the inter atomic 
potential. The main advantage of this method, and the reason for it being used more than the NtT method, is the 
faster convergence of the stress-stress fiuctuation term relative to the convergence of the strain-strain fluctuation term. 
Two steps are required to calculate the elastic constants using this formula. As in the NtT formalism first one has 
to find the zero stress reference mat rix and volume ho and Vq and then run the NVT MD program to calculate the 
elastic constants using equation ( |2.2| ). 



III. VALIDATION AND COMPARISON WITH EXPERIMENT 



In order to validate the programs and to check the potential at various conditions we performed a series of calcu- 
lations of the physical properties of copper using our programs and compared the results with experimental data. 

First, we calculated the thermal expansion up to the melting point at zero external pressure. Our results are 
compared with neutron scattering dataEj in Fig. |l|. As expected, the calculated values are close to the experimental 
values at low temperatures (We must bear in mind that the lattice constant at zero temperature was used as input in 
the fitting of the potential). The agreement remains good even at elevated temperatures, with a maximal difference 
of about 10%. 

Second]#- the phonon spectrum, Lu(k) , was determined by calculating the time Fourier-transform of the following 
quantitiesEZi: 

JV 

fa{k,t) ^"^{ua ■ Vi)cos{k ■ ri) (3.1) 

4=1 

where, fi and Vi are the time dependent positions and velocities of the ith atom. The polarization vectors are defined 
as ui — k/k in the longitudinal case, while Uti and Ut2, the transverse polarizations, are perpendicular to ui and to 
each other. 

We then calculated the phonon spectrum along the (100) and (110) directions. Due to the importance of long-time 
correlations the calculations were made under conditions of constant energy without temperature control. We ran 
the MD program for 20000 time steps with temperature and pressure control, 10000 time steps with temperature 
control only, and then fa{k,t) was saved to disk during 200000 time steps under energy rGcnservation conditions. The 
resulting spectrum, shown in Fig. |[ is in good agreement with neutron scattering datatS. 

Thirdly, to test the algorithms used to obtain the elastic constants, we calculated their variation with temperature 
using both approaches described in section II. The system was initially equilibrated for 50000 steps, and 300000 
subsequent steps were used to calculate the fiuctuations. We found that the convergence of the PR method was 
slower than that of the NVT method. To monitor the convergence of the algorithm, we compared values of the 
elastic coefficients calculated along directions which are symmetrically equivalent. In the final result, after 300000 
time steps, in both methods, the three equivalent elastic constants differ by no more than 10%. From this variance, 
one can estimate the accuracy of the calculation. Comparison with the experimental values measured by neutron 
scatteringtS, shows that the shear moduli of the model-lattice and of real copper soften in the same way within the 
estimated accuracy of the simulations (e.g. 10%). The results are shown in Fig. ||. 

The typical computer time needed to calculate the elastic constants varies between 4 days for a calculation in the 
NVT ensemble to one week in the NtT ensemble on a Pentium-Pro computer. 
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IV. INFLUENCE OF POINT DEFECTS 



After successful validation of the pure system, point defects were introduced into the lattice by either removing an 
atom (vacancy) or adding one (interstitial). The simulations were done using a box having a fixed number of sites. 
Periodic boundary conditions were used throughout. In order to prevent mutual annihilation of the defects, we used 
defects of only one species in each run. We first used simulated annealingl3 to find the most stable configuration 
of a vacancy and of an interstitial in the model. The formation energy of a vacancy was found to be 1.27eV. The 
lowest energy configuration of an interstitial atom was found to be a CipO)-split interstitial with a forniation energy 
of 3.28eV. These formation energies agree well with experimental dataE3, and with previous simulationsEl. Formation 
energies of other possible configurations of interstitials were determined using the conjugate gradient method, and are 
summarized in Table |. We point out that the difference between these energies is on the order of 0.1 eV (^ lOOOK), 
which means that interstitials of various configurations may be present at high temperatures. 

The influence of point defects was determined by calculating the shear elastic coefficients in samples containing 
different concentration of defects using MD. For any concentration of defects, the number of atoms and unit cells in 
the sample is kept constant, so that the number of externally introduced defects is conserved. Under this constraint, 
atoms and defects are free to diffuse, agglomerate etc. We first discuss vacancies. We found that the influence of 
vacancies on the shear elastic coefficients is too small to be resolved within our stated accuracy. For example Fig. 
^ shows that C44 is independent of the concentration of vacancies. This is correct both at low temperatures and 
near the melting point. Previously, influence of vacancies on the elastic properties at zero temperature was treated 
analytically by Dederichs, et alE3 who also reviewed earlier work. The two body potentials used by these authors are 
less accurate than the one used here. They found a reduction of 2-3 % in the values of C44 and 1.5-3 % in values of 
C per percent of vacancies. The effect is indeed small, falling within the uncertainty limits of our calculation. 

We turn now to the influence of interstitials. In contrast to the vacancies, interstitials seem to cause a large 
softening of of the shear moduli, (see Fig. ^). At low temperatures, we can qualitatively compare the results witii 
the analytical calculation of Dederichs et al, made at T=0KEI, and with the irradiation experiment of Holder et a£3, 
made at T=4.2K. In the latter, the defects produced by the irradiation are Frenkel pairs and not interstitials. The 
common trend is that defects induce a large softening of the shear moduli. The softening is anisotropic, meaning that 
C44 softens more than C . Bearing in mind the different conditions under which the different studies were made, the 
qualitative agreement can be regarded as good. 

The main new thrust of this work regards studying the influence of point defects on the shear moduli at increasingly 
higher temperatures, up to the melting temperature Tm- We emphasize that due to the use of periodic boundary 
conditions, the crystal can be super-heated and T^/ is higher than the thermodynamic melting temperatuce. 

In general, our results can be explained by the combined effect of volume change and the dia-elastic effected. Volume 
changes resulting from introduction of defects arise due to the specific volume per defect being larger than that of 
an atom residing at a normal lattice site. We can compare the change of the elastic coefficients caused by defects 
with that induced by heating the crystal up to a temperature where thermal expansion increases the volume by the 
same amount. If the dependence of the elastic coefficients on the volume is the same, it would mean that effect of 
the defects is simply to expand the lattice. We demonstrate such a comparison in Fig. ^ for data obtained at T= 
400K. As one can see in this figure, the change of C due to thermal expansion (without defects) is identical with that 
induced by increasing density of defects at a constant temperature when plotted against the lattice constant. This 
proves that C is affected only by volume changes. On the other hand, the change of C44 due to thermal expansion is 
very different from that induced by increasing the density of defects, and cannot be explained by the volume increase 
alone. We attribute this behavior to the dia-elastic effect. Briefly, the dia-elastic effect describes the coupling of the 
internal vibration modes of the interstitial with the phonons. Some of the internal vibration modes, the so called 
resonance modes, have frequencies within the band of the usual lattice spectrum. Mode coupling will lower the long 
wavelength part of the phonon branch related to the elastic coefficients. For example, the (100) split interstitial is 
coupled to (100) shear strain. If the lattice is deformed in this way the response of the atoms forming the defect 
would be larger thauptiifi-i-esponse of atoms occupying normal lattice sites. The large response leads to a reduction of 
C44 elastic coefficientO'O. On the other hand, C is associated with (110) shear strain, which does not couple to the 
(100) split interstitial. Consequently, the decrease of C' is induced only by volume expansion, as indeed seen in Fig||. 

Our results indicate that at high temperatures we can make one of the shear moduli vanish by adding a sufficient 
number of self-interstitials to the sample. The number of defects needed for this purpose is temperature dependent. 
The absolute change of the shear moduli per % of interstitials is larger for C44 at all temperatures, as summarized in 
Table || 

However, the relative sensitivity of C becomes larger at high temperatures as melting is approached, because its 
initial value at elevated temperatures is already very small. This is plotted in Fig. |[ Consequently,C is the first one 
to vanish. This result is quite surprising, as it is opposite to what one would expect based on the low temperature 
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behaviouiE3. 

We determined the melting temperature within our model at different temperatures and defect concentration, by 
monitoring the structural order parameter after 200000 time steps. The results are summarized in Fig. 0. Each point 
in the figure represents the results of several runs made at constant T and with the same defect density, but at a 
different random initial configuration. The melting process found is mechanical melting, and the sample collapses 
due to the vanishing of C . This was found by correlating the temperatures at which C vanishes with the melting 
temperature. The melting process is rather fast, its duration being on the order of a typical vibration time of the 
crystal. Thus, this is a bulk melting mechanism. As can be seen from our results the interstitials reduce the mechanical 
melting point, and a concentration of about 0.005 lowered the melting point by about 60°K. The same cannot be said 
about vacancies, which did not reduce the mechanical melting point even at very high concentration (up to 4%). 

It turns out that the specific volume per atom at the phase boundary in Fig. |^ is fixed at 13.3lA'^. The same 
specific volume is found by extrapolating the dependence of C on the volume to the value at which C =0 . This 
specific value is consistent with the jSOiecific volume of the liquid phase at T=1356K, the thermodynamic melting 
temperature, as suggested previouslyE3'E3. The way one reaches that specific volume, either by heating or by adding 
defects at a constant temperature does not seem to matter. 

In summary, we have simulated the influence of point defects on the shear moduli of Cu at temperatures up to Tm- 
Mechanical melting occurs due to the fuUfillment of the Born criterion, where one of the shear moduli of the solid 
vanishes. We found that mechanical melting of the crystal is brought about by the vanishing of C , rather than that of 
C44. Melting occurs once the specific volume reaches a critical value. This value is independent of the path by which 
it is reached, either by thermal expansion alone, or by any combination of thermal expansion and of expansion due 
to interstitials introduced into the sample. Since the specific volume per vacancy is smaller than the atomic volume, 
vacancies do not expand the crystal. Hence, vacancies play no role in bulk melting. 
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Configuration 


E^{eV) 


vacancy 


1.27 


(100) split 


3.28 


Octahedral 


3.45 


(110) split 


3.55 


(111) split 


3.78 



TABLE I. Formation energies of various point defects. 



Temp(K) 


AC44{GPa/%inters.) 


AC {GPa/%inters.) 


400 


-39 ± 13 


-4.5 ±2.5 


1400 


-23 ± 1.7 


-3.5 ±2.0 


1450 


-18 ±2.5 


-4.1 ±2.0 


1480 


-15.5 ±1.3 


-7.8 ±3.5 



TABLE II. Change in shear moduli per percent of self-interstitials at various temperatures. 
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FIG. 1. Thermal expansion coefRcient of copper as,function of temperature obtained from our simulations compared with 
data determined by a neutron scattering experimenlEj. 
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FIG. 2. Phonon dispersion crpve at 400K along (100) and (110). Solid lines are the result of the simulations, while the circles 
are experimental results at 80Kll3. 
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FIG. 3. Variation of the shear moduli with temperature: PR method, •; NVT MD, □; experimental dataLj, dashed line. 
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FIG. 4. Dependence of C44 on the concentration of defects at T = 1400^^, showing the difference between the influence of 
vacancies and interstitials. 
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FIG. 5. The variation of the shear moduli vs. the lattice constant. The solid line represents thermal dilation and the symbols 
represent dilation due to interstitials. The dashed line is a guide to the eye. 
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FIG. 6. Relative softening of the shear moduU per percent of interstitials as function of temperature. The Unes are guides 
to the eye. 



1580 



1530 



1480 



Melt 



Solid 



1430 

0.002 0.004 0.006 

Concentration of split-interstitials 

FIG. 7. The influence of self interstitials on the melting temperature obtained using periodic boundary conditions (a crystal 
without free surfaces). 
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